#!/usr/bin/perl
# 2011/12/26 1000bpのwindowのスコアを比較するためのデータを作成
# 2012/01/05 1kb non-overlapping windowでscoreを比較する

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);

use File::Basename;
$| = 0;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [wig file #a] -b [wig file #b] -o [output image]\n";
}

my %opt;
getopts("a:b:o:", \%opt);
my $WIG_a = $opt{a};
my $WIG_b = $opt{b};
my $out_image = $opt{o};
(my $TITLE_a = basename($WIG_a)) =~ s/\.wig//;
(my $TITLE_b = basename($WIG_b))=~ s/\.wig//; ;

my %map_a = &readMap($WIG_a);
my %map_b = &readMap($WIG_b);
sub readMap{
	my ($file) = @_;
	my %data;
	my $fh;
	if($WIG_a =~ /\.gz$/){
		$fh = IO::File->new("gzip -dc $file |") or die "cannot open $file: $!";
	}else{
		$fh = IO::File->new($file) or die "cannot open $file: $!";
	}
	my $chr = '';
	while($_ = $fh->getline()){
		s/\r?\n//;
		if(m/^track/){
			next;
		}
		if(m/chrom=(chr\d)/){
			$chr = $1;
			next;
		}
		my ($pos, $value) = split /\t/;
		my $cate = int($pos/5000) * 5000 + 500;

		if($value < 0){
			next;
		}

		$data{"$chr\t$cate"}{'sum'} += $value;
		$data{"$chr\t$cate"}{'num'}++;
	}
	$fh->close();
	return %data;
}

my $tmp = 'tmp.txt';
my $fh_out = IO::File->new($tmp, 'w') or die "cannot write $tmp: $!";
$fh_out->print("a\tb\n");
foreach my $key(keys %map_b){
	unless(exists $map_a{$key}{'sum'}){
		next;
	}
	$fh_out->printf("%f\t%f\n", $map_a{$key}{'sum'} / $map_a{$key}{'num'}, $map_b{$key}{'sum'} / $map_b{$key}{'num'});
}
$fh_out->close();


my $Rcommand = <<END;
R --vanilla --slave <<EOF
	data <- read.table("$tmp", header=T, sep="\t");
	attach(data)


	pearson = cor(a,b,method="pearson")
	kendall = cor(a,b,method="kendall")
	spearman = cor(a,b,method="spearman")

	jpeg(file="$out_image", width=400, height=400)
#	pdf(file="$out_image")
	plot(b, a, pch=20, cex = .8, xlab="$TITLE_a", ylab="$TITLE_b")
	abline(lm(b ~ a), col="red")
	legend("topleft", pch = 20, legend=c(paste("pearson = ", sprintf("%.3f", pearson), sep=""),
	paste("kendall = ", sprintf("%.3f", kendall), sep=""),paste("spearman = ", sprintf("%.3f", spearman), sep="")))
	graphics.off()
#	dev.off()

	cat(paste("pearson = ", sprintf("%.3f", pearson), ", kendall = ", sprintf("%.3f", kendall), ", spearman = ", sprintf("%.3f", spearman)), sep="")
EOF
END

my $cor = `$Rcommand`;
print "$cor\n";
unlink $tmp;



